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ABSTRACT 

We present self-consistent high-resolution simulations of NGC4038/4039 (the "Antennae galaxies") 
including star formation, supernova feedback and magnetic fields performed with the JV-body/SPH 
code Gadget, in which magnetohydrodynamics are followed with the SPH method. We vary the 
initial magnetic field in the progenitor disks from 10~ 9 to 10~ 4 G. At the time of the best match 
with the central region of the Antennae system the magnetic field has been amplified by compression 
and shear flows to an equilibrium field value of ~ 10 /xG, independent of the initial seed field. These 
simulations are a proof of the principle that galaxy mergers are efficient drivers for the cosmic evolution 
of magnetic fields. We present a detailed analysis of the magnetic field structure in the central 
overlap region. Simulated radio and polarization maps are in good morphological and quantitative 
agreement with the observations. In particular, the two cores with the highest synchrotron intensity 
and ridges of regular magnetic fields between the cores and at the root of the southern tidal arm 
develop naturally in our simulations. This indicates that the simulations are capable of realistically 
following the evolution of the magnetic fields in a highly non-linear environment. We also discuss 
the relevance of the amplification effect for present day magnetic fields in the context of hierarchical 
structure formation. 

Subject headings: methods: JV-body simulations — galaxies: spiral — galaxies: evolution — galaxies: 
magnetic fields — galaxies: kinematics and dynamics 



1. INTRODUCTION 

Within the framework of the hierarchical galaxy for- 
mation picture, galaxies assemble and ev olve via merg- 
ers of smaller progenitor g alaxies (e.g. iWhite fe Reesl 
[19781 IWhite fc Frenld[l99lT i. Thus, galaxy interactions 
are essential for the understanding of structure forma- 
tion. In the bottom-up-picture of structure formation 
dwarf galaxies merge to form the first galaxies at an early 
epoch of the universe. Later, there is still a continuous 
merging of fully evolved galaxies. The further growth 
of galaxies progresses through a combination of mergers 
and diffuse accretion of gas. 

Interactions of g alaxies change their dynamics 
drast ic ally (see e.g. iToomre fe Toomrd 119721 IBarnesI 



1992. iHernquist fe Barnes! 
Burkert fe Naabl 120031 " 



1994 IBarnesI 11999 and 
Gonzalez-Garcia et al.l 120061 ) as 
the gravitational potential is changing rapidly during 
the interaction. Since the gas component is dissipative 
and most sensitive to changes of the gravitational 
potential, it is strongly affected during the interaction 
and driven to the galaxy centers, eventually causing 
bursts of star formation (iBarnes fe Hernauistl I 1992 L 
Mihos fe Hernauistl jl99j. iBarnes fe Hernauistl Il996i 



Bekki fc Shibval UM&L ISpringell 12000. IBarnesI [2002 



Cox ct al. 2006, Naab ct al. 200 



Bournau d et al.l|2005l.v .... 

Robertson et all 120061 iCox et al.l 120081 iHopkins etakl 
2008). So far simulations of interactions and mergers 



of disk galaxies have only been investigated with re- 
spect to changes in stellar dynamics, gas flows, star 
formation (SF ) or formation o f cent ral supermassive 
black holes (|Di Matteo et al.l 120051 ISpringel etakl 

kotarba@usm.lmu.de 



I 2005U ISpringel et al.l [2005a[ IRobertson et af] 
2006L iDeBuhr et alj 120091 I Joh ansson et al.l T2009b. 
I Johansson et al.ll2009af h However, the dramatic impact 
of mergers on the gas flows will directly affect the 
magnetic fields of the systems (and vice versa) via the 
induction equation of magnetohydrodynamics (MHD) 
and the Lorentz force. The magnetic fields will change 
their morphology following the motion of the gas and 
will be amplified by shocks and gas inflow. 

Changes in the magnetic field structure, on the other 
hand, might influence gas flows, local collapse and the 
morphology as well as the star formation activity. Lo- 
cal, interacting galaxies are the perfect laboratories for 
investigating the effects associated with their magnetic 
fields. However, the timescales for galaxy mergers are 
far too long to observe these processes directly. The only 
observational possibility to study the time evolution of 
mergers is to consider different systems at different evo- 
lutionary stages. However, the available sample of in- 
teracting nearby galaxies is too small to investigate the 
evolution of magnetic fields in detail. Thus, numerical 
simulations pose a promising tool to study the magnetic 
field evolution in interacting systems. 

The structure of an interacting system strongly de- 
pends on the properties of the progenitor galaxies. Thus, 
matching observed nearby interacting systems with sim- 
ulations in space and time can give us an idea of the 
properties of their progenitors, e.g. their sizes, gas frac- 
tions and relative velocities. Furthermore, comparing 
simulated systems with observations helps to asses the 
performance of the applied numerical method. Numer- 
ical methods supported by these comparisons can then 
be used to study processes in the early universe, when 
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galaxies were very different from present-day galaxies. 
High resolution simulations of the formation of indi- 
vidual galaxies i n a full cosmological context (see e.g. 
iNaab et alll2007l ) including magnetic fields could help us 
in understanding the processes leading to the magneti- 
zation of the universe. This type of study would comple- 
ment earlier semi-analytical studies that investigated the 
possibility of magnetic field seeding by gala ctic winds in 
a cosmological context (jBertone et al |[2006T) . 

The standard theory of magnetic field amplification 
in galaxies is the so called Galactic Dynam o based 
on the mean field theory (see iKulsrudl 119991 for a re- 
view) . Within this theory, turbulent motions of the ion- 
ized gas driven by stellar activity lead to the genera- 
tion of a random magnetic field (a-effect). This ran- 
dom magnetic field (particularly its radial component) 
can then be amplified by the differential rotation of 
the galaxy (O-effect), leading to an efficient dynamo 
action which results in an exponential growth of the 
magnetic field (the afi-dynamo). However, dynamos 
may probably only work efficiently if magnetic helic- 
ity is transported away from the differe nt ially rotating 
disc ([Brandenburg fe Subramani an 2005). Grcssc l et al.l 
(2008) performed high-resolution box-simulations which 
demonstrate that a dynamo may operate if supernova 
explosions release magnetic helicity from the disc. How- 
ever, for an efficient magnetic helicity transport out 
from a galactic disk, galactic winds or galactic foun- 
tains may be required. This might be a problem par- 
ticularly for massive galaxies due to the deeper poten- 
tial well. The fact that it is difficult to get an efficient 
dynamo is generally addressed as the dynamo problem. 
Different sol utions, e . g. tu rbulence driven by large-scale 
SN-bubbles (lFerrierelll992T ) or the Cosmic Ray Dynamo 
(jHanasz et al.1 120091 ) have been proposed. These solu- 
tions describe the exponential growth of a small-scale 
magnetic seed field which is amplified up to present-day 
values within several Gyr. However, recent observations 
indicate that magnetic fields in galaxies have been al- 
ready very strong (comparable to present-day galactic 
magnetic fields) at very high redshifts, at a time when 
the u niverse was only t ss 6 Gyr old (zwl) (jBernet et al.1 
20081). Former obs ervations of damped Ly-a systems by 
Wolfe et all ([1992) indicate that progenitors of galactic 
discs had magnetic fields of a few /iG even at z ps 2 (t s» 3 
Gyr). The very fast amplification required to generate 
the strong magnetic fields at high redshifts can probably 
not be achieved with any Galactic Dynamo model (see 
e.g. lArshakian et al.l [20091 ) . Thus, alternative possibil- 
ities for the amplification of galactic ma gnetic fields on 
shorte r timescales need to be explored. iLesch fc Chibal 
(|1995t ) have shown analytically that strong magnetic 
fields in high redshift objects can be explained by the 
combined action of an evolving protogalactic fluctuation 
and el ectrodynamic p roces ses providing magnetic seed 
fields. iWang fc Abell ()2009f ) performed numerical simu- 
lations of the formation of disc galaxies within an col- 
lapsing halo imposing a uniform initial magnetic field of 
f 0~ 9 G. The initial field was amplified by three orders of 
magnitude within approximately 500 Myr of evolution. 
The amplification might be due to the combined effects 
of magnetic field compression during the collapse and 
amplification of the unifor m initial field by d ifferential 
rotation as studied also in iKotarba et al.l (|2009l ). These 



studies indicate, that the amplification of magnetic fields 
might be a natural part of the galaxy formation process. 
However, interactions of galaxies, which were more fre- 
quent at earlier times, pose another promising possibility. 

Although it would be worthwile to consider cosmo- 
logical studies of structure formation including mag- 
netic fields in the long run, numerical studies of in- 
teracting magnetized systems in the local universe 
may serve as a first step towards a more complete 
scenario. These studies help us in understanding 
the morphological e volution of galaxies, their star 
forma t ion histories dSpringel et al.l l2005d ICox et all 
20081 iBournaud et all 120071 iDi Matteo et al.l I2008L 
Jesseit et al.ll2009l INaab fc Ostrikerll2009( ). and as we will 
show in this paper also their magnetic histories. The sys- 
tem NGC 4038/39, also known as the Antennae galaxies, 
is one of the best examples for an ongoing local merger. 
It is also the by far best observed interacting galaxy sys- 
tem. 

In this paper we present further steps towards a more 
complete numerical representation of the Antennae sys- 
tem as a prototype for interacting galaxies. For the first 
time we will follow the evolution of the magnetic field 
in a galaxy interaction simulation. We also address the 
general question whether smoothed particle magnetohy- 
drodynamics (SPMHD) is capable of following the evolu- 
tion of magnetic fields in interacting systems at reason- 
able accuracy. In a previous paper we have shown that 
SPMHD is well suited for following the evolution of mag- 
netic fields in isolated disk galaxies (jKotarba et al.1120091 ) 
so the study of interacting systems is a natural extension 
of this earlier study. 

The paper is organized as follows: Section [2] summa- 
rizes the properties of the Antennae system as known 
from observations and theory. In section [3] we describe 
our numerical methods fsection l3.1[) . the setup of the iso- 
lated disks (section I3.2.1[) and the match to the observed 
Antennae system (section |3.2.2[) . A detailed analysis of 
the evolution of the system is presented in section 13.31 
where we discuss the evolution of the magnetic field (sec- 
tion l3.3.1~T) the numerical stability of our simulations (sec- 
tion I3.3.2j) and the self-regulation of the magnetic field 
amplification (section 13.3 . 3[) . In section 2] we describe 
our method to calculate artificial radio maps (section 
14. ip and present applications to the isolated disk and 
the Antennae simulations ( section [4. 2p . The artificial ra- 
dio maps derived from the simulations can be compared 
against the radio observations of the system, thus provid- 
ing a further tool for constraining the numerical model 
and method. We conclude in section [5] and briefly dis- 
cuss the relevance of our simulations in the context of 
hierarchical structure formation. 

2. PROPERTIES OF THE ANTENNAE SYSTEMS 

The Antennae system is relatively nearby, the es- 
timated distances range from 13 to 25 Mpc. The 
smaller distances are favored by m ethods based on pho - 
tometry of the red giant branch ([Saviane 
whereas the larger distances are estimated from ob- 
servations of the Type la superno va 2007sr in the 
southern tail ( Schwciz er et al.l 120081 ). But note also 
that sometimes even values up to d = 29 Mpc have 
been adopted in the li terature ([Fabbiano et al.l I2001L 
iZezas fc Fabbianol I2002D . In this paper, we apply the 
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Fig. 1. — Total synchrotron emission (contours) and magnetic 
field vectors of polarized intensity at 4.86 GHz based on VLA data 
(yellow), overlaid on a DSS image (blue - white background) (Dig- 
itized Sky Survey, Palomar and UK Schmidt telescopes). The con- 
tour levels are 0.005, 0.12, 0.30, 0.53, 1.2, 2.1, 3.3, 5.3, 9.0, 17 and 
24 mj y /beam-area. The resolution is 17"xl4". Credit: IChvzvl 
(2005) 



conventional distance of 22 Mpc for all relevant flux cal- 
culations. Given the large variety of high quality ob- 
servations (iWhitmore et al.lll999t iNeff fe UlvestadlfepOC' 



Wilsc 


n et al. 2000, Hibbard et al. 2001, Chvzv & Beck 


2004, 


Wang et al. 2004, 


Brandl et al.l 20051 iZezas et al. 


2006, 


Brandl et all 2009 


) several authors tried to find 



initial condi tions for simulations representing the Anten- 
nae system. iToomre fc Toomrel (|1972h first presented re- 
stricted three-body simulations which already explained 
the formation of tidal arms and bridges as a result of 
tidal interaction during the merger. Follow-up investi- 
gations confirmed the early results by studying the de- 
tailed galactic dynamics using self-con sistent, multiple- 
component galaxy models (Barnes 1988). Further studies 
added star-formation (jMihos et al.lll993l ) to the modeling 
process and aimed at constraining the influence of dark 
halo mass profile s on the development and morphology 
of the tidal tails (iDubmski et al.lll996h . 

Recently, iKarl et al.l (|2010| ) developed a new model of 
the system, not only focussing on its plane-of-sky ap- 
pearance, but also on fitti ng the observed line-of-sight 
velocity structure (see also IKarl et al.ll2008l). This study, 
alongside with new observations (jZhang et al.l2010f) , sug- 
gests that the localized intense starburst sites observed in 
the overlap-region can be explained as the imprint of the 
interpenetrating process of the two merging disks follow- 
ing their second encounter. These results contrast with 
previous numerical simulations which find the current or- 
bital phase of the Antennae system to lie somewhere be- 
tween the first and the s econd c l osest encounters (e.g. 
Toomre fc Toomrd 119721 iBarnesI 119881 iDubinski etaLl 
199G). 

In this paper we use the model of lKarl et al.l (|2010| ) and 
focus on the central region of the Antennae system and 
its magnetic fields. Fig. [TJ shows a DSS image (Digitized 
Sky Survey (Blue) , conducted with the Palomar and UK 
Schmidt telescopes and digitized by the Catalogs and 
Surveys Group of the Space Telescope Science Institute) 



of the central region of the Antennae system, i.e. the 
galactic disks and bases of the tidal tails. Overlaid are 
contours of total radio synchrotron emission (tracing the 
total magnetic field energy) and magnetic field vectors 
(derived from polarized intensity). The strength of the 
magnetic field is 20 fiG on average. The highest values of 
more than 30 (j,G are reached in the overlapping regio n 
and the centers of the galaxies (|Chvzv fc Beck! 120041) . 
Thus, the magnetic field is roughly twice as strong as 
the typical values of 5 to 10 obs e rved in isol a ted spi- 
ral g alaxies (e.g. IBeck et al.l 119961 , IBeckl [20071 IKrausd 
2009). On the upper left (east), there is a large region 
with highly ordered magnetic field lines, most probably 
tracing the gas flow at the root of the lower (southern) 
tidal tail. This gas is also visible as a HI ridge which 
exten ds far out along the southern tail (Hibbard et al. 
120011 ). The structure of the magnetic field most likely 
traces the recent gas motions induced by tidal forces dur- 
ing the merger. Apparently, not much of the magnetic 
field structure of the progenitor galaxies have survived 
the interaction. The progenitors were presumably typi- 
cal spirals with a spiral magnetic field pattern (see e.g. 
[BeckT[2009l ). As the magnetic field is tightly linked to 
the motion of the gas, the structure of the field in a sys- 
tem which has recently undergone a violent interaction 
should mainly resemble the recent kinematic evolution. 
It does not depend on long-term processes like the Galac- 
tic Dynamo, which i s believed to be important in isolated 
spiral galaxies (e.g . IBra ndenburg & Subram anianl 1 2005 



Gressel et al. l 12008. Beck 2009, Dubois fc Tevssierl I2009L 
Elstner et al.l I2009l iGiss inger et al. 20091 lHanasz et all 
20091 ). In other words, nonlinear systems lose the mem- 
ory of their initial conditions. Hence, numerical studies 
of the kinematics of merging systems including magnetic 
fields should be able to represent the observed magnetic 
fields in nearby interacting systems also without includ- 
ing long-term processes. 

3. SIMULATIONS 

3.1. Numerical methods 

All simulation s were perform ed with the N-body/SPH- 
code Gadget (|Springell 120051) . Gravitational interac- 
tions between the par ticles are evaluated w ith a hier- 
archical tree method ()Barnes fc Hull H986I) . The dy- 
namics of Lagrangian fluid elements are followed using 
a SPH formulation which conserve s both energy and 
entropy ([Springel fc Her nquist 2002) including the evo- 
lution of magnetic fields wh i ch wa s implemented and 
tested bv iDolag fc Stasvszvnl ([2009). The code has al- 
ready been used to investigate the evolution of ma gnetic 
fields in isolated spiral galaxies (jKotarba et al .1120091 ) and 
to compare different implementations of the SPH for- 
mulations and implementations in the SP H Code VINE 
(jWetzstein et al J 120091 fNelson et all 120091) . These stud- 
ies demonstrated the importance of a sensible treatment 
of the numerical divergence of the magnetic field (V • B) 
in SPH simulations, as it can lead to artificial magnetic 
field growth. An implementation utilizing Euler poten- 
tials, which by constructio n poses a V • B-free prescrip- 
tion of magnetic fields (see lPrice fc Bate! |2007) for more 
details) circumvents this problem. However, using the 
Euler potentials, the magnetic field is essentially mapped 
on the initial particle arrangement. Thus, if the initial 
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configuration significantly changes shape during the sim- 
ulation, regions carrying conflicting values of Euler po- 
tentials (i.e. values, which do no longer allow for a fi- 
nite and unambiguous calculation of their gradients) can 
overlap and the ability of the Euler potentials to repre- 
sent the magnetic field correctly is lost. This can lead to 
defective magnetic field calculations, especially in kine- 
mat ically vigorous syste ms like interacting galaxies (see 
also iBrandenburgj [20091 . Therefore, all simulations pre- 
sented in this paper have been performed using the stan- 
da rd (direct) mag netic field implementation. In contrast 
to iKotarba eTalT l|2009h we now also apply the Lorentz 
force and artificial magnetic dissipation applying an ar- 
tificial magnetic dissipation constant of as — 0.5. The 
latter does not only allow for magnetic field redistribu- 
tion and reconnection, but also lowers the numerical di- 
vergence as it helps to smooth artificially high magnetic 
fields arising from intrinsic constraints of the numeri- 
cal prescription. In this sense, it poses a regularization 
scheme similar to smoothing of the magnetic field. Both 
schemes reduce th e numerical noise and V • B errors by 
a similar amount (Dolag & Stasyszyn 2009). However, 
the dissipation scheme is based on physical considera- 
tions, whereas the smoothing scheme is completely ar- 
tificial. Thus, using the dissipation scheme, the factor 
hV ■ B/|B| (where h is the so called smoothing length, 
which poses the typical length scale of spatial derivatives 
in SPH calculations) is restricted to a value of approx- 
imately unity. Values of order unity have been shown 
to be low enough to guarantee a physically meaningful 
evolution of the magnetic fields in SPH simulations, par- 
ticularly preventing artificial magnetic field growth. This 
threshold is actually defined by simulations using Euler 
potentials, for which the numerical divergence measure 
hV • B/|B| is of order unity alt hough the physical di ver- 
gence i s zero by definition (see IKotarba et all 120091 and 
section 15.3.21 for more details). 

Further more, we do no t use a viscosity limiter as sug- 
gested by iBalsaral (|1998l ). because applying this limiter 
resulted in an increased growth of the magnetic field. 
This is most likely a numerical artefact, as the limiter 
lowers the viscosity in regions of strong shear flows, thus 
suppressing velocity diffusion and leading to a higher ve- 
locity dispersion and higher velocity gradients, which in 
turn lead to artificia lly enhanced magnetic field growth 
(|Kotarba et al.H2009t ). 

All simulations are performed including radiative cool- 
ing assuming a primordial gas composition together with 
a homogeneo us and time-independen t extragalactic UV 
background (|Haardt fc Madaul flggG) . We include star 
formation and the associated supernova feedback, but 
exclude explicit supernova-driven galactic winds, fol- 
low ing the sub-resolution multi phase model developed 
by ISpringel fe Hernquistl (120031). in which the ISM is 
treated as a two-phase medium ([McKee fc Ostrikerlll977l 
Uohansson fc Efstathioull2"006[ ): Cold clouds are embed- 
ded in a tenous hot gas at pressure equilibrium. Stars 
form from the cold clouds in regions were n > n t h = 
0.128 cm~ 3 with the shortlived stars supplying an energy 
of 10 51 ergs to the surrounding gas by supernovae. The 
threshold density, n t h, is determined self-consistently in 
the model by requiring that the equation-of-state (EOS) 
is continuous at the onset of star formation. The param- 
eters governing the model (see Tab. [T]) are set to produce 



TABLE 1 
Parameters of initial setup 



Disk parameters 


total mass 


-Mtot 


1.34 x IO^Mq 


disk mass 


M disk 


0.075 Mtot 


bulge mass 


-Mbulgc 


0.025 Mtot 


mass of the gas disk 


Mgas 


0.2 M disk 


exponential disk scale length 


Id 


8.44 kpc 


scale height of the disk 


I'D 


0.2 l D 


bulge scale length 


Id 


0.2 l D 


spin parameter 


A 


0.1 


virial velocity of the halo 


«vir 


160 km s" 1 


half mass radius 


^half 


«12 kpc 


half mass circular velocity 


"half 


Ri249 km s" 1 


half mass rotation period 


Thalf 


Ri295 Myr 


initial magnetic field 


B 


10~ 9 to 10" 6 G 


Multi-Phase model parameters 


gas consumption timescale 


tMP 


8.4 Gyr 


mass fraction of massive stars 


Pmp 


0.1 


evaporation parameter 


A 


4000 


effective SN temperature 


Tsn 


4 x 10 8 K 


cold cloud temperature 


Tec 


1000 K 



a star formation rate of ~ 1M© yr" 1 for a Milky Way- like 
galaxy in isolation. 
The implement ation used in this p aper has been 



tested in detail (ISpringel et all 120011 ISpringell 120051 
ISpringel et aTll2005bl . iDolag fc Stasvszvnl 120091 ) and ful- 



fills the esta blished requirements for nu merical methods. 
Particularly, IDolag fc Stasvszvnl (|2009h have shown that 
the MHD-implcmentation performs well in various test 
problems, incl uding different shock t ube problems, the 
Fast Rotator (IBalsara fc Spicerlll999f) . the Strong Blast 
(e.g. [Balsara fc Spicerlll999D and the Orszag-Tang Vor- 
tex (|Orszag fc Tangjll979fh 

3.2. Setup 
3.2.1. Isolated disks 

The Antennae system has most likely formed through 
the interaction of two formerly isolated spiral galaxies. 
In this section we present the properties of the isolated 
progenitor model disks used in our simulations. The ini- 
tial conditions for the spi ral galaxies are realize d using 
the method described bv | Springel et all (|2005b[ ) which 
is based o n iHernquistl (| 19931 ). The galaxies consist of a 
iHernquisH (| 19901) profile cold dark matter halo, a rota- 
tionally supported expone ntial stellar disk, an exponen- 
tial gas dis k and a stellar | Hernauistl (|1990[ ) bulge com- 
ponent (see lKarl et al.ll2010l for more details). The halo, 
stellar disk and bulge particles are collisionless N-body 
particles. The gas is represented by SPH particles. 

Possible initial conditions resulting in a good fit of 
the present-day properties of the Antennae galaxies have 
been tested in a large parameter study by iKarl et all 
(2010) (see also section 15. 2.2j) . In this paper we present 
results using the initial condition parameters of this 
study which result in the best match to the central region 
of the Antennae system. The parameters describing the 
initial conditions of the two galaxies can be found in Ta- 
ble [TJ Particle numbers and softening lengths are listed 
in Table [2j The disks are kinematically s table if evolved 
in isol ation as has been shown in detail in IKotarba et all 
( 2009). In the following we thus only address the evolu- 
tion of the magnetic fields. 

For simplicity, the initial seed magnetic field is assumed 
to be homogeneous with only one non-vanishing compo- 
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TABLE 2 

Particle numbers and softening lengths 



Component initial particle number fixed gravitational 

softening length e [pc] a,b 



Halo 4.0 X 10 5 80/h 

Disk 4.8 x 10 5 20/ft 

Bulge 2.0 X 10 5 20/h 

Gas 1.2 X 10 s 20/h 

Stars 20/h 

Total 1.2 x 10 6 



Fig. 2. — Gas surface density S at time t = Myr (upper panel) 
and t = 400 Myr (lower panel) , overlaid with magnetic field vectors 
for the simulation with Bq = 10 -6 G. The length I of the vectors is 
normalized to a minimal value B m i n = Bo / \/2 and displayed loga- 
rithmically according to I = 5 log (B/B m i n ), i.e. I = corresponds 
to B B m in or smaller and I = 5 to B = 10-B m ; n . 



nent of B x = Bq. This choice is justified, as it takes 
more than one Gyr of dynamic evolution until the present 
plane-of-sky-appearance of the system has developed in 
our simulations. Thus, the particular structure of the ini- 
tial magnetic field should not be of significance for the fi- 
nal result. We use two different values for the initial field, 
B = 10~ 9 G and Bo = 10 ~ 6 G for the isolated galaxies, 
and additionally two intermediate values, Bo — 10~ 8 G 
and Bo = 10 G for the Antennae simulations. The 
smallest value of Bo — 10 -9 G is the typical value of 
the observed intergalactic magnetic field £?igm (see e.g. 
iKronberg et~aLll2008H and the highest, B a = 10~ 6 G, is 
motivated by the typical value of several /LtG observed in 
spiral galaxies. As much larger or much smaller values 
are not observed, these values cover the range of real- 
istic initial fields. However, we have also performed a 
simulation of the Antennae system with an initial mag- 
netic field value of B = 10~ 4 G in order to study the 
physical behaviour of the system in an extreme situation. 
We do not include neither large-scale dynamo processes, 
nor turbulent motions on scales smaller than w 100 pc 
which are not resolved in our simulations. The mean 
velocity dispersion a = \J (lP)~Zn. — (v)n. n . (where the 
mean is taken over the nearest 64 ± 5 neighbors within 
the smoothing kernel) during the isolated disc simula- 
tion is of order of 5 km s" 1 with approximately 30% of 
the particles having dispersions > 5 km s _1 and only 
a few percent > 10 km s _1 . These values are some- 
what lo wer than the values f ound in recent grid simula- 
tion s bvlWang fc Abell ll2009h ( see Fig. 3 of their paper) 
and lAgertz et al.l (|2009f ). They both find typical disper- 
sion values of approximately 10 km s _1 in their com- 
parable disc galaxy simulations. However, these authors 
use lower temperature floors for the dense gas component 
found in the star-forming regions, resulting in a clumpier 
disc structure and thus probably in an enhanced turbu- 
lence in the hot diffuse component of their discs. This 
may explain the discrepancy in the measured velocity dis- 
persions. Since the dispersion values in our simulations 
are rather low, we do not expect any significant amplifica- 
tion of the magnetic field in the isolated galaxies. Conse- 
quently, the magnetic field gets only redistributed during 
the simulation, developing a spiral pattern as the differ- 
ential rotation continues to wind it up (Fig. [2]), while the 

overall value of \B\ remains of order Bo throughout the 
simulation. After the magnetic field has been wound up 
by differential rotation, it is highly ordered in the disc 
region (r > 5 kpc) and more " turbulent " in the inner re- 
gion of the galaxy. Thus, the inner magnetic field is not 
visible in the lower panel of Fig. [2] due to the averaging 
calculation of our plotting routines. 

The evolution of the absolute value of the magnetic 
field as a function of time is shown in Fig. [3] for Bo = 
10~ 6 G (red line) and Bo = 10~ 9 G (orange line). In the 
beginning of the simulation, the initially homogeneous 
magnetic field gets wound up and thereby amplified due 
to the differential rotation and ass ociated shear flows b y 
roughly a factor of two (see also iKotarba et al.l [2009) . 
After approximately one half mass rotation period the 
magnetic field has been redistributed to a mostly toroidal 
pattern by the winding process, and the amplification 
ceases. In the subsequent evolution, the strength of the 
field decreases slowly due to magnetic dissipation, which 



a The minimum SPH smoothing length for the gas particles is l.Oe. 
The Hubble constant is assumed to be /i - 0.71 in this paper. 
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Fig. 3. — Btot = y + By + B% as a function of time for the 
Antennae simulations with an initial field of Bo = 10 -9 G (green 
line), B = 10 -8 G (blue line), B = W' 7 G (dark blue line), 
B = 10" 6 G (black line) and B = 10" 4 G (black dashed line), 
respectively, and for the progenitor disks simulations with Bo = 
10 — 9 G (orange line) and Bo = 10 — 6 G (red line), respectively. The 
magnetic field of the isolated disks does not evolve significantly. 
For the mergers the field is amplified to lOfiG independent of 
the initial field strength in the disks. 

is the only process causing magnetic field diminution in 
our simulations. The velocity dispersion of the gas also 
leads to the development of a z-component of the mag- 
netic field (not shown) which, however, remains smaller 
than all other components by more than one order of 
magnitude throughout the simulation. In summary, the 
magnetic field gets redistributed to form a spiral pattern 
(Fig [2]) and retains on average its initial value through- 
out one Gyr of evolution. This behaviour is the same for 
both the weaker and the stronger initial magnetic field. 

Fig. [4] shows the mean numerical divergence /iV-B/|B| 
as a function of time in isolated galaxy simulations (red 
line) with Bq = 10~ 6 G. The mean was taken over three 
simulations with the magnetic field in the plane of the 
disk and inclined as in the setup of the Antennae sim- 
ulation (see section I3.2.2|) . respectively. Although the 
numerical noise increases with time, it remains clearly 
below the tolerance value of unity (see also section l3.3.2p . 

The SFR in the isolated disks is roughly constant 
throughout the simulations (not shown). Starting at a 
value of approximately 2 Mq yr _1 and then decreasing 
slightly to approximately 1.7 Mq yr _1 after 1.3 Gyr of 
evolution due to gas consumption. There is no significant 
difference in the evolution of the SFR compared to the 
same simulation without any magnetic field, indicating 
that the presence of a global magnetic field of order 10 -6 
G or lower does not affect the gas flow enough to hinder 
or abet the collapse of gas. 

3.2.2. The match to the Antennae system 

The simulations presented here arc taken out of a suit 
of self-consistent simulations designed as a large param- 
eter study to match the morp hological and kin ematical 
properties of the Antennae (see lKarl et al.ll2~010D . In this 
study, we initially set two equal-mass galaxies, each mod- 



FlG. 4. — /iV ■ B/|B| as a function of time for the Antennae 
simulations with an initial field of Bo = 10 -9 G (green line), Bo = 
10" 8 G (blue line), B = 10" 7 G (dark blue line), B = 10~ 6 G 
(black line) and Bo = 10~ 4 G (black dashed line), respectively, and 
mean divergence for isolated simulations with Bo = 10 — 6 G (red 
line) . The values stay below the tolerance value of unity throughout 
the simulation in every run. 

eled as in section [3.2.11 and residing in its own dark mat- 
ter halo, on nearly-parabolic Keplerian two-body orbits 
with given ellipticity e, pericenter distance r p , and initial 
separation r scp . The disk orientation in the orbital plane 
is given by a pair of angles (t, w), which, for each galaxy, 
specify the adopted inclination with r espect to the orbital 
plane and the pericentric argument (|Toomre fc Toomrd 
[1971 ). There is no hot gas component surrounding the 
galaxies initially. The initial field is assumed to be ho- 
mogeneous with only one non-vanishing component of 
B x = Bq. After the simulation has finished we deter- 
mine the time of best match, the viewing direction of 
the observer, a common center-of-mass, and a distance 
scale factor C in order to create a mock observation which 
can be compared to p rojections of the HI data cube from 
Hib bard et all (|200lh . If the result does not prove satis- 
factory up to a level admissible by optical inspection, the 
simulation is repeated choosing a different set of initial 
parameters. Several key parameters regarding the ellip- 
tical orbit, the relative orientation of the galaxy disks, 
and the internal structure of the progenitor galaxies are 
varied in order to find the best match (for details, see 
iKarl et aT1l2010l ). The final parameters used in this study 
are shown in Tab. [3l 

Starting on their initially set orbit, both galaxies evolve 
corresponding to their isolated evolution (section I3.2.1[) 
until they reach the point of their first closest approach 
it s» 650 Myr). At this time, the prominent tidal arms, 
which we use as tracers for the dynamical history of 
the encounter, start to develop. On the other hand, 
the detailed structure of the galactic main bodies can 
only be seen in our simulations resulting from the recent 
splash during the second encounter (t 1180 Myr). The 
time of best fit, i.e. the time, at which the simulation 
matches the appearance of the Antennae system in the 
sky and the observed line-of-sight velocities, is reached 
at t B M « 1250 Myr. 
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TABLE 3 

Antennae simulation parameters 



Initial orbit parameters 


disk orientation 


NGC 4038 


NGC 4039 


i 


60° 


60° 




30° 


60° 


ciiipticity 


e 


0.96 


pericenter distance 


r P 


7 kpch- 1 


initial separation 


fscp — ^vir 


160 kpch -1 


Analysis parameters 


time oi best match 




1.25 Gyr 


direction to observer a 


(0,^,0) 


(208,282,0) 


distance scale 


c 


2.0 



a The viewing direction is specified by a series of ro- 
tations in the following order around the x-, y-, and 
z-axis. 



3.3. Evolution of the Antennae system 
3.3.1. Magnetic field evolution 

We have run several simulations using the setup de- 
scribed in section 13.2.21 The initial field was again as- 
sumed to be homogeneous with only one non-vanishing 
component of B x = Bq at the beginning of the simula- 
tion. We performed simulations with five different initial 
magnetic field strengths of B = 1CT 9 , 10~ 8 , 1(T 7 1CT 6 
and 10 -4 G for comparison. 

Fig. [5] shows the line-of-sight magnetic pressure 
-fmag = B 2 /8tt in the simulation with Bq = 10~ 6 
G at different time steps, overlaid with contours of 
the stellar surface density £ s tars- The particle data 
has been transferred to a grid of 80x80 cells using 
the TSC procedure (Triang ular Shaped Cloud, see e.g. 
iHocknev fc Eastwood! Il988f i. As it takes roughly 600 
Myr (i.e. approximately two half mass rotation periods) 
before the first encounter of the galaxies, the magnetic 
field has enough time to redistribute and form a realistic 
configuration in each of the galaxies prior to the merger 
(upper left and right panel, see also section 13.2.11 and 
Fig. [2]) . The formation of the tidal arms is visible in the 
stellar density distribution but also in the distribution 
of magnetic pressure (central right and lower left panel). 
At time of best fit (lower right panel) most of the gas 
has been driven into the central region of the Antennae 
system. Thus, the magnetic pressure reaches its highest 
values in this region. 

The temporal evolution of the absolute values of the 
magnetic fields for the simulations with different initial 
field values is shown in Fig. [3] (black dashed, black, dark 
blue, blue and green lines). In all cases (except for the 
run with Bq = 10~ 4 G), similarly to the simulations of 
the isolated galaxies (red and yellow lines), we see a mild 
amplification of the initial magnetic field in the beginning 
of the simulation due to the winding process. However, 
as the disks are not oriented parallel to the xy-plane, this 
initial amplification is slightly weaker than in the isolated 
disks. The reason is that the initial magnetic field now 
does not lie in the plane of the disks and thus the radial 
component of the magnetic field is weaker compared to 
the simulations of the isolated galaxy. In the case with 
the weakest initial field the magnetic field gets ampli- 
fied by more than two orders of magnitude during the 
interaction, whereby the most violent amplification oc- 
curs during the first encounter at t sa 650 Myr. In the 
case with Bq — 10~ 6 G, however, the amplification is 



relatively modest. The evolution of the magnetic field 
for the simulation with the highest initial field (dashed 
line) is different: At the very beginning of the simula- 
tion, the high magnetic overpressure drives the gas out 
of the galaxies, thus "blowing" them up. Consequently, 
the magnetic field decreases by one order of magnitude 
within 100 Myrs due to attenuation and continues to 
decrease until the first encounter. At the time of the 
first encounter, it is only very weakly amplified. At the 
time of best match, the value of the magnetic field is 
approximately 10 /iG within the numerical precision, in- 
dependent of the initial seed field. This is roughly half 
the value derived from observations. The origin of this 
discrepancy might be observational as well as numerical 
and will be briefly addressed in section 14.21 

3.3.2. Numerical stability 

Fig. |4] shows the arithmetic mean of the numerical di- 
vergence /iV • B/|B| as a function of time for the Anten- 
nae simulations with the different initial magnetic field 
strengths (black dashed, black, dark blue, blue and green 
lines) . For each simulation, there is an increase of the di- 
vergence during the first encounter, whereby the value of 
the numerical divergence increases with decreasing ini- 
tial magnetic field. This is not surprising: If the mag- 
netic tension is strong enough to overcome the gas pres- 
sure, the Lorentz force acts on the particles in a way that 
magnetic tension is released. On the other hand, if the 
magnetic pressure is significantly weaker than the gas 
pressure, chaotic motions of the particles driven by the 
encounter can fold the magnetic field on small scales - as 
small as the smoothing length - until the magnetic ten- 
sion becomes dominant. This leads to a more irregular 
magnetic field and a higher numerical divergence. Thus, 
the numerical divergence is lowered in the presence of a 
stronger magnetic field. 

Fig. [5] shows the mean line-of-sight numerical diver- 
gence /iV • B/|B| in the simulation with Bq — 10~ 6 G 
at different time steps, overlaid with contours of the stel- 
lar surface density £ s tars to indicate the morphology of 
the galaxies. The particle data has been transferred to 
a grid using the TSC procedure as before in Fig. [5j 
Before the first encounter (central left panel), there are 
regions of high (of the order of 10) numerical divergence 
at the "edges" of the galaxies (upper panels, compare 
also Fig. [5]). This high numerical divergence measures 
can be ascribed to defective SPH calculations in these 
regions. The particle density there decreases to zero due 
to the vacuum boundary conditions (which are usually 
used in this type of simulations). Thus, the particle 
distribution within one smoothing length changes rather 
abruptly. Some SPH operators, including the divergence 
operator, are not well sampled in such a situation, lead- 
ing to high numerical errors in theses estimators. As 
soon as the particle distribution is smoothed out as a 
consequence of the interaction (central left to last panel, 
compare also Fig. [5]), this effect vanishes. However, it 
is only a small fraction of particles which are affected by 
this defective calculation. Thus, the arithmetic mean of 
the numerical divergence is lower in the beginning of the 
simulation than after the first encounter (Fig. @|. Fur- 
thermore, comparing Fig. [5] with Fig. [5] shows that in 
regions with the highest magnetic field values the numer- 
ical divergence is relatively low. 
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Fig. 5. — The Antennae simulation with an initial field of Bq = 10 6 G. Colors visualize the line-of-sight magnetic pressure .Pmag 

= B 2 /8tt 

(in units of 10 — 18 g cm -1 s — 2 ) and contours correspond to stellar surface density E stars . The contour levels are 0.005, 0.02, 0.08, 0.32, 
1.28, 5.12 and 20.48 M Q pc" 2 . 
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Fig. 6. — The Antennae simulation with an initial field of Bq = 10 6 G. Colors visualize the mean line-of-sight numerical divergence 
/iV • B/|B| and contours correspond to stellar surface density £ a t a rs- The contour levels are 0.005, 0.02, 0.08, 0.32, 1.28, 5.12 and 20.48 
M Q pc" 2 . 
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We have performed the same simulation with an ini- 
tial magnetic field of 10~ 6 G but without applying the 
Lorentz force (not shown). In this simulation the mag- 
netic field got amplified extremely violently by orders of 
magnitude to clearly unphysical values after the first en- 
counter. The magnetic field got amplified much above 
the maximal value seen in the simulations presented 
above, and did not converge. This behaviour shows, that 
it is actually the Lorentz force, i.e. the backreaction of 
the magnetic field on the gas, which constrains the am- 
plification. The unrealistic violence of the amplification 
can be traced back to the high V • B values of several 
hundreds developing in this simulation. However, as ap- 
plying the Lorentz force helps to lower the divergence in 
SPMHD simulations, these results are not surprising. 

This can also be seen in simulations including the 
Lorentz force, but starting with an initially very weak 
magnetic field. We have performed an additional simula- 
tion with Bq = 10 -20 G (not shown). In this simulation, 
the divergence grew to a maximal value of 2.5 during the 
first encounter, subsequently dropping again to values 
below unity. The magnetic field - and thus the Lorenz 
force - was very weak in this simulation, nevertheless, the 
divergence was still lowered to values of order unity. The 
magnetic field was amplified by ten orders of magnitude 
to a value of 10 -10 G during the first encounter, which 
is still four orders of magnitude lower than the maximal 
value seen in Fig. [31 This demonstrates that one can not 
start with an arbitrary low magnetic field and end up 
at micro-gauss levels after the first encounter. However, 
as the subsequent interaction between the two galaxies 
drives further turbulence, the magnetic field continued to 
grow after the first encounter at a rate of approximately 
one order of magnitude per 100 Myr. Thus, at a time 
of 1.1 Gyr, the magnetic field reached a value of w 10~ 6 
G and slowly converged towards the maximal value seen 
in Fig. |3J Note that during this steady growth of the 
magnetic field the numerical divergence was actually de- 
creasing. 

In the simulation with Bq = 10~ 4 G, the numerical di- 
vergence measure grows up to a value of approximately 
0.4 already at the beginning of the simulation. This is be- 
cause the high magnetic pressure "blows up" the galaxies 
and thus excites strong turbulent motions which in turn 
result in a more irregular magnetic field. We note that 
this behaviour shows that the value of the numerical di- 
vergence mainly depends on the irregularity of the mag- 
netic field, which is also the reason why the numerical 
divergence remains relatively small in the simulations of 
the quiescent evolution of the isolated galaxies (red line). 
This can also be understood theoretically: The numerical 
SPH divergence-operator calculates the weighted sum of 
the differences of the magnetic field of a particle and its 
neighbouring particles within a smoothing length. Thus, 
the higher the degree of irregularity of the magnetic field, 
the higher the numerical divergence. However, the nu- 
merical divergence should not be mistaken for a physical 
divergence, as it is only a measure of numerical small- 
scale (i.e. smaller than one smoothing length) fluctua- 
tions of the field. This can be seen in simulations using 
Euler-Potentials, where the physical divergence is zero 
by definition, but the numerical divergenc e has still val- 
ues of order unity (see lKotarba et al1l2009( ) for a tangled 
magnetic field. Thus, lowering the numerical divergence 
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Fig. 7. — /iV ■ B/|B| as a function of the total magnetic field for 
the Antennae simulation with Bq = 10~ 9 G at time of the first 
encounter (t fti 650 Myr). Grey dots correspond to the values of 
each particle, the solid line is the mean value for a given magnetic 
field strength. The values of the numerical divergence are widely 
distributed over the range of magnetic field strength. 
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Fig. 8. — Same as Fig. but at the time of best fit (f-BM ~ 
1250 Myr). The values of the numerical divergence are widely 
distributed over the range of magnetic field strength and even lower 
for higher magnetic field values. 



below this tolerance value of unity should be sufficient 
to guarantee a physically meaningful evolution of the 
magnetic field. Using the Euler-Potentials in the Anten- 
nae simulations would most probably result in a much 
weaker amplification of the magnetic field which should 
not be considered physical, because Euler-Potentials are 
not suitable for simulations of kinematically vigorous sys- 
tem (see section ETTj) . 

Moreover, the applied SPH implementation is geared 
to ensure that the numerical divergence measure does not 
alter the evolution equat ions for the magnetic field (see 
iDolag fc Stasvszvnl[200l . Thus, even if the divergence 
operator measures a numerical divergence, it does not 
influence the m agnetic field evoluti o n dire ctly. This has 
been shown by lPrice~fc Monaghan (2005), who demon- 
strated that a magnetic monopole can be advected with- 
out causing numerical instabilities. It can also be seen 
comparing Fig. [3J and Fig 0J The lower the initial mag- 
netic field, the weaker the magnetic field shortly after 
the first encounter, although the numerical divergence is 
higher for lower initial fields. Thus, there is no direct 
dependance of the magnetic field strength on the value 
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of the numerical divergence. Fig. [7] and [5] show the nu- 
merical divergence as a function of the magnetic field 
strength at the time of the first encounter and at the 
time of the best fit, respectively. We show the plots for 
the simulation with Bq = 10~ 9 G, as the amplification of 
the magnetic field is the most efficient and the numerical 
divergence is the highest in this simulation. Thus, a pos- 
sible dependance of the magnetic field on the divergence 
measure should be the best visible in this simulation. 
However, there is no significant correlation, and the val- 
ues of the numerical divergence are widely distributed 
over the range of magnetic field strength. At the time 
of best fit, they are even lower for higher magnetic field 
values (see also Fig. [6]). This behaviour is qualitatively 
the same for all initial magnetic field values. Of course, 
the amplification is more efficient for lower initial fields, 
thus one could argue that it is the amplification efficiency 
which depends on the numerical divergence value. How- 
ever, in the beginning of the simulation with Bo = 10~ 4 
G the magnetic field actually decreases with increasing 
numerical divergence showing that non- vanishing numer- 
ical divergence not necessarily leads to an amplification of 
the magnetic field. Rather, the amplification efficiency 
is restricted by the strength of the Lorentz force: The 
higher the magnetic field, the stronger the Lorentz force 
braking the motions which lead to an amplification of the 
magnetic field. Thus, the lower the initial field, the more 
efficient its amplification. Hence we conclude that as long 
as the numerical divergence remains as low as the numer- 
ical divergence seen in simulations with Euler-Potentials 
(i.e. lower than unity), the evolution of magnetic fields 
in SPH simulations is physically meaningful. 



3.3.3. Self -regulation of the amplification 

The magnetic field is expected to get enhanced through 
field line compression in shocks and field line stretching 
in shear flows. However, in the framework of MHD, any 
motion of gas leading to an amplification of the magnetic 
field will be suppressed by the magnetic field itself via the 
Lorentz force as soon as the magnetic energy gets com- 
parable to the kinetic energy of the gas. The magnetic 
energy is then converted into kinetic energy of the gas, 
thus maintaining equipartition between the magnetic and 
gas kinetic energy density, or equivalently, the magnetic 
and the hydrodynamic pressure -Phyd = l/2pv 2 . In par- 
ticular, the magnetic field is expected to be in equi parti- 
tion with the turbulent ener gy of the gas (see e.g. iBeckl 
120071 and IChvzv et ;il.ll2l)l)~i. as only velocity gradients 
can lead to an amplification of the magnetic field via 
the induction equation. Thus, the self-regulation of the 
strength of the magnetic field seen in our simulations can 
be ascribed to equipartition between the turbulent and 
magnetic pressures. In order to analyze this behavior, we 
have examined the central region of the system in more 
detail, and also performed a comparison simulation with- 
out magnetic fields. We define the turbulent pressure as 
-fturb = 1/2/CfturbJ with the turbulent velocity of the i-th 
particle defined as 



Uturb(i) = 




where k = x,y, z and 

i v k (j) — v k (i) 

<b(o = y 2 ^ =1 % — (2) 

with TV = 64 ± 5 being the number of the nearest neigh- 
bors. 

We briefly note that the thermal pressure does not di- 
rectly affect the evolution of the magnetic fields. Accord- 
ing to the induction equation of MHD, the magnetic field 
evolution is determined by the velocity field alone. 

Fig. [9] shows from left to right and top to bottom 
the gas number density n, the turbulent pressure Pturb, 
the stellar surface density S s tars, the velocity dispersion 
a (calculated as before in section 13.2. ip . the magnetic 
pressure P mag and /3 = Pturb/Pmag in the inner region 
(innermost 28 kpc) of the system at time of best fit for 
the simulation with Bq = 10~ 6 G. This value of Bq 
is comparable to the typical magnetic field value ob- 
served in spiral galaxies, which is why we have chosen 
this simulation for our analysis. Assuming a distance of 
22 Mpc this region w ould comprise approximately 4.37'. 
IChvzv k. Beck1l2004l use a distance of 19.2 Mpc and ob- 
serve an area of « 3.5', scaled to a distance of 22 Mpc this 
gives an area of approximately 22.4 kpc across. Thus, our 
model has a bigger extent by a factor « 28/22.4 = 1.25, 
which still is in qualitative agreement with observations 
(see section [2]) . The particle data has been transferred 
to a spatial grid using the TSC procedure and aver- 
aged over the z-direction (i.e. the line-of-sight) with 
z € [—14 kpc, 14 kpc] and z — defined as the center 
of mass of the system. The turbulent and magnetic pres- 
sures are given in units of 10~ 15 g cm -1 s~ 2 , correspond- 
ing to 6.242 x 10 -4 eV cm -3 , i.e. the highest values are 
approximately 100 eV cm~ 3 . 

In order to be able to recognize whether the magnetic 
field itself has a significant effect on the turbulent pres- 
sure in the system, we have applied the same analysis to 
a simulation without magnetic fields. Fig. [10] displays 
from left to right the mean line-of-sight gas number den- 
sity n, turbulent pressure Pturb and stellar surface den- 
sity S s tars in the inner region (innermost 28 kpc) of the 
not magnetized system at time of best fit (calculated as 
before in Fig. [9]). 

Comparing Fig. [9] and Fig. 1 101 shows that in the simu- 
lation with magnetic fields (Fig. |9]) the gas distribution 
is more compact, whereas the turbulent pressure distri- 
bution is "disrupted". Particulary in the northern (up- 
per) galaxy the turbulent pressure distribution is more 
extended in the magnetized case (Fig. [9]) than in the 
simulation without magnetic fields (Fig. ITU]) . These dif- 
ferences probably develop because the gas is more likely 
to move along magnetic field lines than perpendicular to 
them and thus the velocity distribution is altered. The 
gas distribution is in both cases more extended than the 
stellar distribution, with the stellar density being highest 
in the centers of the galaxies (« 10 3 M Q pc -2 ). The stel- 
lar distribution is not significantly changed in the pres- 
ence of a magnetic field. Since the galaxies have a low 
gas fraction (20%), the total gravitational potential is 
dominated by the stellar component in the inner region 
of the Antennae system. Thus, the distribution of the 
high density gas (> 10 cm~ 3 ) is almost unaffected by 
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Fig. 9. — The Antennae simulation with an initial magnetic field of Bo = 10~ 6 G. From left to right and top to bottom: Mean line-of-sight 
gas number density n, turbulent pressure Pturb = V^P^turb' s *- euar surface density Estars, gas velocity dispersion <r, magnetic pressure 
Pmag = B 2 /8n and f) = Pt U rb/Pmag in the inner region (innermost 28 kpc) of the system at time of best fit (tr3M ~ 1250 Myr). 



log( n [10 5 cm 3 ] 



l°g( Pt^b [1CT 15 g sT 2 ] ) 



log( S sl . rs [M s pc- 2 ] ) 




x [kpc] 



Fig. 10. — The Antennae simulation without including magnetic fields. From left to right: Mean line-of-sight gas number density n, 
turbulent pressure Pturb = V^P^turb ant ^ s t eiiar surface density £ s tars m t ne inner region (innermost 28 kpc) of the system at time of best 
fit (*bm ~ 1250 Myr). The turbulent pressure is highest in the overlapping region between the two merging galaxies. 
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Fig. 11. — Temporal evolution of Pturb (black line), P ny d (blue 
line), P,y (green line) and P ma g (red line) for gas particles with a 
number density > 0.005 cm -3 in the simulation with Bq = 10 — 6 G. 
Pturb an d Pmag are of a comparable order of magnitude throughout 
the simulation and almost in equipartition at time of best fit. 

the presence of the field. As stars are expected to form 
in high density regions, it is not surprising that the dis- 
tribution of stars formed in our simulations is also inde- 
pendent on the presence of the field. In the magnetized 
case (Fig. [HI, the gas velocity dispersion a (lower left 
panel) is of the order of 10 to 20 km s _1 within the 
galaxies. The distribution of magnetic pressure (lower 
central panel in Fig. [S]) is slightly different compared to 
the distribution of the turbulent pressure (upper central 
panel): The highest turbulent pressures are reached in 
the centers of the galaxies, whereas the magnetic pres- 
sure is highest in the overlapping region of the galaxies. 
Moreover, there is a ridge of magnetic pressure at the 
root of the southern tidal tail (indicated by the black 
oval) which is not visible in the distribution of turbulent 
pressure. This differences most probably originate in the 
magnetic field being a vector instead of a scalar quantity. 
A fully random magnetic field is not amplified efficiently 
by isotropic compression. Thus, only in regions with 
strong shear flows which stretch and therefore straighten 
the magnetic field it can be amplified efficiently. How- 
ever, the energy range of the magnetic pressure is overall 
comparable to the energy range of the turbulent pressure. 
Thus, j8 = Pturb/fmag (lower right panel) is in the rage 
1 to 10 almost everywhere, which means that the mag- 
netic pressure is of the order of the turbulent pressure or 
slightly lower. 

Fig. fTTI shows the temporal evolution of the turbulent 
pressure Pturb ~ u turb (black line), the hydrodynamic 
pressure Phyd ~ v 2 (blue line) the magnetic pressure 
Pmag ~ B 2 (red line) and the "dispersion pressure", cor- 
responding to the velocity dispersion, i.e. P a — l/2pa 2 
(green line) for gas particles with a number density 
> 0.005 cm -3 in the simulation with B = 10~ 6 G. The 
hydrodynamic pressure is higher than the magnetic pres- 
sure by roughly three orders of magnitude throughout 
the simulation, which should be expected from theory as 
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Fig. 12. — Same as Fig. 1111 but for the simulation with Bq = 
10 -4 G. Pmag is much higher than Pt ur b m the beginning of the 
simulation, but decreases to the level of equipartition within 400 
Myr. Pt ur b and P ma g are almost in equipartition at time of best 
fit. 

it is not the value of the velocity itself, but the veloc- 
ity gradients which determine the evolution of magnetic 
fields. The turbulent, dispersion and the magnetic pres- 
sures are of the same order of magnitude until the first 
encounter (except of the beginning of the simulation). 
After the encounter, the turbulent and the dispersion 
pressures are always slightly higher (by a factor of ap- 
proximately five and ten, respectively) than the magnetic 
pressure. At time of best fit the turbulent and magnetic 
pressures are again of the same order of magnitude, as 
already indicated in the last panel in Fig. [S] The fact, 
that the magnetic pressure never exceeds the turbulent 
pressure indicates, that the magnetic field amplification 
is restricted to magnetic field values corresponding to 
the equipartition level between turbulent and magnetic 
pressure. This is exactly what is expected from theory 
and explains the self-regulated saturation of the magnetic 
field strength in our simulations (Fig. [3J 

Fig. [12] shows the same quantities as in Fig [TT] but 
for the more extreme simulation with Bq = 10~ 4 G. In 
the beginning of this simulation, the magnetic pressure is 
three orders of magnitude higher than the turbulent pres- 
sure (because this initial magnetic field is two orders of 
magnitude higher than the expected equipartition value 
of several /xG and P ma g oc B 2 ). Within the first 50 Myr 
of evolution the magnetic pressure drops by one order of 
magnitude. Simultaneously, the turbulent and dispersion 
pressures increase by the same amount. This is because 
the high magnetic pressure "blows up" the galaxies in 
the very beginning of the simulation and thus drives a 
lot of turbulent (or chaotic) motions. After the first 50 
Myr, the difference between the turbulent and the mag- 
netic pressure is only one order of magnitude and the 
system is able to relax again. Thus, the turbulent and 
dispersion pressures start to decrease, and the magnetic 
pressure continues to decrease further. After approxi- 
mately 400 Myr the magnetic pressure is of the order of 
the turbulent and dispersion pressures. Shortly before 
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the first encounter, it has reached a value slightly be- 
low the turbulent pressure. In the subsequent evolution, 
similar to the simulation with Bq = 10~ 6 G, the mag- 
netic pressure always stays below the turbulent pressure. 
However, the evolution of the pressure components is al- 
tered compared to the simulation with Bq = 10 ~ 6 G. 
Particularly, the second encounter (visible as a tempo- 
rary increase of the pressure values) preceding the time 
of best fit in Fig. [TT]is shifted by approximately 100 Myr 
to later times in the simulation with Bo — 10~ 4 G. This 
difference develops because the strong magnetic field in 
the beginning of the simulation with Bq = 10 ~ 4 G alters 
the gas distribution significantly and thus changes the 
evolution of the whole system. In summary, this com- 
parison clearly shows that interacting galactic systems 
always tend to reach equipartition, independent of the 
initial ratio of magnetic to turbulent pressure. 

As already discussed in section l3.3.2l without applying 
the Lorentz force the magnetic field gets amplified much 
above the value of equipartition between magnetic and 
turbulent pressure, and does not converge. Thus, it is 
actually the Lorentz force, i.e. the backreaction of the 
magnetic field on the gas, which yields the self-regulation. 

Finally, we compared the SF rates in the simulations 
with different initial magnetic field strengths with the 
SF rate in a simulation without any magnetic field (not 
shown). The SFR after the first encounter in the simu- 
lation with Bq = 10~ 6 G showed to be slightly lower (by 
a factor of approximately two) than in the simulation 
without or with a weak magnetic field, indicating that 
the presence of the magnetic field hinders the collapse of 
gas. However, this influence is not strong enough to alter 
the SF history significantly. 

4. SIMULATED RADIO EMISSION AND POLARIZATION 

MAPS 

4.1. Computation method 

In order to compare our results directly with observa- 
tions, we compute artificial radio emission and polariza- 
tion maps from our simulation data. For this purpose, 
the magnetic field components and the stellar density 
from the SPMHD simulations have been again trans- 
ferred to a three-dimensional grid. The following calcu- 
lat ions have bee n performed with an IDL code developed 
bv iWiatrl (|2006ft . The calculations of the total and po- 
larized synchrotron intensity and the calculation of the 
polarisation angle have been performed in the standard 
way a cco rding the to the following fo rmulae fsee lLongairl 
119971 and lRvbicki fc Lightmanl l7986 for more details): 

The total synchrotron emission J v Eit a given frequency 
v is given by 



J v = 



4 a T c 3 



3 /i 



' 27rm 3 



2e 3 



iVB^+nA- 



(3) 



where the magnetic field B is the only input from our 
simulations. The frequency v and the index of the power 
spectrum of the relativistic cosmic ray (CR) electrons 
p are input parameters. The latter is assumed to be 
2.6 in this paper, c orresponding to the value given by 
iChvzv fc"l3 cck 200^. The constants are the Thompson 



cross section ot — 0.665 x 10~ 24 cm 2 , the magnetic per- 
meability no = 1 (in CGS units), the speed of light c 
and the electron mass m e « 9.1 x 10~ 28 g. The constant 
normalization factor k of the cosmic ray energy spectrum 
can be derived for a given total CR energy Ecu via 



E C r = k 



E 1 ~ p dE. 



(4) 



The observed CR energies in the Milky Way follow a 
steep spectrum from 10 9 to 10 20 eV, whereby super- 
nova remnants (SNR) are the most likely source for 
CRs with energies < 10 18 eV. CRs with higher en- 
ergies may be produced in Jets of pulsars or black 
holes, and are probably of extragalactic o rigin (see 
iThe Pierre AUGER Collaboration et al.l 12008ft . Given 
the steep fall-off of CR abundance with energy we assume 
an energy range of E min = 10 9 eV to E max = 10 15 eV 
in our calculations. Furthermore, as CRs in this energy 
range are most likely produced in SNRs, we assume the 
CR distribution to be proportional to the stellar density, 
with a typical value of the mean specific energy density 
of eq B. = 1 MeV m~ 3 for CR protons (see e.g. iFerrierel 
12001ft . However, we apply a cutoff at an energy den- 
sity of ecR = 100 MeV m~ a . The energy density of CR 
electrons is roughly 100 times lower than the energy den- 
sity of CR protons, thus the mean energy density for CR 
electrons is assumed to be 10 keV m~ 3 . 

J„ is calculated within every grid cell at a frequency 
of v = 4.86 x 10 9 Hz (corresponding to the observed 
frequency). The total intensity /tot of the synchrotron 
radiation is subsequently obtained by integration of the 
emission along the line-of-sight. 

The degree of polarization n of any electromagnetic 
radiation is defined as the amount of its polarized inten- 
sity I p compared to the amount of its total intensity Jt t- 
The synchrotron emission of a single radiating charge 
is always polarized elliptically, because the light compo- 
nent for which polarization is parallel to the magnetic 
field projected onto the plane of sky (In) has a differ- 
ent refraction index than the perpendicular component 
(/]_). However, as charges gyrate along the magnetic field 
lines, the elliptical components will cancel, as emission 
cones will contribute equally from both sides of the line- 
of-sight. Thus, for any reasonable distribution of parti- 
cles that varies smoothly with pitch angle, the radiation 
will be partially linearly polarized and thus characterized 
by the terms 7m and I± . The degree of linear polarization 
for particles of a single energy can then be expressed as 



1%) 



Ix{v) + h(yY 



(5) 



where Itot{v) = I±_{v) + 1\\ (v)- If the energy spectrum of 
the radiating particles follows a power-law (here N{E) = 
nE~ p dE), it can be shown that 



n = 



p 



i 



p- 



(6) 



Thus, in a homogeneous magnetic field, the degree of 
polarization is very high (approximately 73% for p — 
2.6). However, when integrated along the line-of-sight, 
opposite polarization cancels out and the observed degree 
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of polarization is therefore usually much lower than the 
theoretically expected value. 

The polarized intensity I p depends on the Stokes pa- 
rameters Q and U according to 



with 



(7) 

(8) 
(9) 



where the integration is performed along the line-of-sight 
and ijj is the polarization angle, defined as the angle be- 
tween the electric field vector of the radiation perpen- 
dicular to the magnetic field (E±) and the a:- axis in the 
xy-plane (the plane of sky), i.e.: 



i P = VQ 2 + u 2 , 

Q = U J J u cos (2tp)ds, 
U = Tl J v sin (2ip)ds, 
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Finally, the observed degree of polarization is 
TT - Ip 

itot 



(10) 

(11) 

(12) 



All calculated values are subsequently convolved with 
a telescope beam c orresponding to the 17 " x 14" beam in 
the observations of IChvzv fc Beckl (|2004l ). i.e. the beam- 
diameter is approximately 1.5 kpc at the distance of the 
Antennae system (assuming a distance of 22 Mpc). The 
shape and sensitivity of the beam is specified by a 2D 
gaussian function. 

4.2. Applications 

Fig. [T3] shows an simulated face-on radio map of 
the isolated disk at t = 400 Myr for the simulation 
with Bq — 10~ 6 G. The particle data of a domain 
with x € [-20 kpc, 20 kpc], y G [-20 kpc, 20 kpc] and 
z G [—10 kpc, 10 kpc] (with the zero-point defined as the 
center of mass of the system) has been transferred to a 
spatial grid with 60 x 60 x 30 cells. Thus, the displayed 
domain comprises 40 x 40 kpc 2 . The colours correspond 
to the stellar surface density, overlaid with contours of to- 
tal synchrotron power. Magnetic field lines derived from 
calculations of polarization are shown in black. To ac- 
count for the spatial isotropy of the emission from any 
emitting volume element, the total flux has been multi- 
plied by the factor 



fobs 



beam 



AlT-d 2 



(13) 



with d = 22 Mpc the distance to the observer and rbeam 
the assumed radius of the beam. Thus, the artificial flux, 
given in mJy, corresponds to what is expected to reach 
the earth from the distance of the Antennae system. 

As already discussed above, the initially homogeneous 
magnetic field gets redistributed by the differential rota- 
tion of the disk, thus developing a spiral pattern which is 
clearly visible in the total emission. The magnetic field 




Fig. 13. — Face on view of the isolated galaxy. Colors correspond 
to the logarithm of stellar surface density (in units of 10 — 5 Mq 
pc — 2 ), overlaid with contours of total synchrotron power. The 
contour levels are 0.001 to 0.01 mJy in ten equally spaced steps. 
Magnetic field lines derived from calculations of polarization are 
shown in black. 



lines trace this spiral pattern. Altogether, the structure 
of the magnetic field is similar to what is observed in typ- 
ical disk galaxies. A similar result has been also obtained 
independently by Kulcsza-Zvdzik ct al. (2009), who have 
performed 3D MHD simulations of barred spiral galaxies 
using a grid code. 

Interestingly, the distribution of the magnetic field 
lines derived from the polarization calculations (Fig. II 3[) 
does not extend as far out in the disk as the magnetic 
field itself (Fig. |2J|. This difference occurs because we can 
only observe polarization where enough CR particles are 
present, which is not the case in the outer parts of the 
galaxy. However, the structure of the magnetic field is 
comparable. 

Fig. [14] shows a simulated radio map of the inner region 
of the Antennae system for the simulation with Bq = 
10~ 6 G at the time of best match (t BM ~ 1-25 Gyr). The 
particle data of a domain with x G [—14 kpc, 14 kpc], 
y G [—14 kpc, 14 kpc] and z G [—56 kpc, 56 kpc] (with 
the zero-point again defined as the center of mass of the 
system) has been transferred to a spatial grid with 75 x 
75 x 300 cells. The total flux was again corrected to the 
isotropy of the emission by multiplying by the factor / bs- 
The contour levels of total synchrotron emission have 
been chosen to be as given in Fig. 3 in IChvzv fc Beckl 
(2004). Thus, they are the same levels as displayed in 

Fig. m 

Given the fact that our simulations are fully self- 
consistent, the similarity between the simulated and the 
observed system is astonishing. The spatial extent and 
distribution of the total synchrotron flux compares very 
favorably with the observations. Also, the highest values 
of total synchrotron emission are reached in the overlap- 
ping region and at the centers of the interacting galax- 
ies. Furthermore, two ridges of ordered magnetic field 
lines, one reaching from one galaxy to the other along 
the overlap region, and one corresponding to the root of 
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Fig. 14. — Inner region (innermost 28 kpc) of the simulated An- 
tennae system. Colors correspond to the stellar surface density 
(in units of 10 — 5 Mq pc — 2 ), overlaid with contours of total syn- 
chrotron power. The contour levels are 0.005, 0.12, 0.30, 0.53, 1.2, 
2.1, 3.3, 5.3, 9.0, 17 and 24 mjy. Magnetic field lines derived from 
calculations of polarization are shown in black. The simulated sys- 
tems compares very well to the observed system (Fig. [!}■ 



the southern tidal tail, naturally develop in our simula- 
tion. However, there are also several differences: There 
is a lack of magnetic fields in the southern tidal tail, i.e. 
the ordered magnetic field structure is not as prominent 
as in the observations, which may be caused by a lack of 
CRs in this region. Also, there is no western spiral arm in 
NGC4038 (the upper galaxy) , which is probably because 
the spiral structure of the progenitors in our simulation 
is not pronounced enough. Furthermore, there is too 
little polarized emission in the outskirts of the galaxies 
and the overlap region is shifted north (down) compared 
to the observations. Moreover, the pitch angle of the 
magnetic field in the isolated galaxy (Fig. [T3|) is rather 
small. The latter can probably be explained by the ab- 
sence of a dynamo process in this simulation. Despite 
these differences the satisfactory match between obser- 
vation and our simulation is encouraging. Thus, our nu- 
merical method already seems to capture the most essen- 
tial basic processes relevant in investigating interactions 
of magnetized galaxies. In particular, our model of the 
Antennae system seems to provide a fair description of 
how this system may have formed. 

A further discrepancy between observation and our 
simulation is the value of the magnetic field itself. The 
magnetic field strength in our simulations saturates at a 
mean value of roughly 10 /iG, only 10% of the simulated 
particles carry magnetic field values of \B\ > 20^G and 
only 1% have \B\ > 50/iG. On the other hand, the mean 
magnetic field strength derived from observations of syn- 
chrotron radiation is approximately 20 /iG. However, the 
observed value is derived assuming equipartition between 
the CR energy and the energy of the magnetic field, an 
assumption which does not necessarily have to hold. Fur- 
thermore, the assumed CR energy density and the mag- 
netic dissipation factor in our calculations are only ap- 
proximate estimates. Given these uncertainties, a differ- 



ence by a factor of two between the observed magnetic 
field and the field strength in our simulations is admissi- 
ble. 

5. CONCLUSION AND DISCUSSION 

We have presented the first fully self-consistent JV- 
body/SPH simulations of the interacting Antennae 
galaxy system including magnetic fields. We show that 
weak magnetic seed fields in the isolated disk galaxies are 
amplified by the gravitational interaction throughout the 
two galactic encounters. Thereby the magnetic pressure 
saturates at a level corresponding to equipartition be- 
tween the turbulent and the magnetic pressure, indepen- 
dent of the initial field strength. Particularly, magnetic 
fields with an initial value higher than the equipartition 
value diminish during the evolution, demonstrating that 
the state of equipartition is the natural state for mag- 
netized galactic systems. An analysis of artificial total 
synchrotron emission and polarization maps provides a 
convincing agreement with the observations. Summariz- 
ing, the method of iV-body/SPH simulations including 
magnetic fields reproduces quite conclusively the compli- 
cated dynamics of the amplification and spatial design of 
magnetic fields in interacting galaxies. 

Moreover, a detailed discussion of the numerical di- 
vergence of the magnetic field in SPH simulations has 
been presented in section 13.3.21 Our analysis strongly 
suggests that numerical divergence measures which are 
smaller than a certain threshold can be considered as 
measures of sub-resolution fluctuations which do not af- 
fect the overall evolution of the magnetic field. Consid- 
ering simulations using the Euler-Potentials, which pose 
a V • B-frec prescription by definition, thi s threshold can 
be as sessed to be /iV • B/|B| « 1 (see also lKotarba et all 
I20M) . 

What can we learn from these simulations for the 
global evolution of cosmic magnetic fields? Within the 
framework of standard CDM hierarchical clustering mod- 
els the formation of large disk galaxies as well as elliptical 
galaxies is characterized by more or less intense merging 
of smaller galactic subunits, e.g. dwarf galaxies, collaps- 
ing gas clouds or globular clusters. If we assume that 
at least some of the accreted subunits have been magne- 
tized by stellar activity (e.g. supernova explosions, stel- 
lar winds or T-Tauri-jets), the merging of such subunits 
to larger galaxies must have been accompanied by a sig- 
nificant amplification and restructuring of the magnetic 
field on galactic scales. The amplification and ordering 
of small-scale magnetic fields to a toroidal configuration 
during th e evolution of isolat ed gal axies was recently 
shown by Ha nasz et al.l (12009D and IDubois fc Tevssierl 
(2009) independently. lHanasz et al.1 (|2009( ) considered 
an axially symmetric galactic disk in which stellar seed 
fields were amplifi e d by a cosmic ray driven dynamo. 
IDubois fc Tevssieil (|2009f) demonstrated the amplifica- 
tion and ordering of small-scale fields seeded by SF activ- 
ity in the context of the formation of a dwarf galaxy with 
significant galactic winds. Complementary to these find- 
ings, our simulations prove that amplification via non- 
axisymmetric three dimensional gravitational interaction 
alone may provide an alternative channel for galactic as 
well as intergalactic magnetic field evolution. In other 
words, given that the structure formation is character- 
ized by a galactic bottom-up architecture, we would ex- 
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pect that within one or two Giga- years the Universe has 
been globally magnetized by the combination of dynamo 
action in isolated galaxies and dynamical amplification 
by interacting galactic objects. However, dynamo ac- 
tion is supposed not to be very efficient in dwarf galax- 
ies since their differential rotation is not strong enough 
(|Gressel et al.l l2008). Thus, at an early epoch of the uni- 
verse, when most of the galaxy population consists of 
dwarfs, magnetic field amplification due to interactions 
may be even more significant. 

With their study of the formation of dwarf galax- 
ies including ma g netic fields and galactic winds, 
iDubois &: Tevssieil (|2009T ) demonstrate an alte rnative 
scenario based on the ideas of lBertone et al.l (|2()()(ii which 
they call the " Cosmic Dynamo" . According to their find- 
ings, galactic winds from young dwarf galaxies eject mag- 
netic field energy into the intergalactic medium, lead- 
ing to a mean intergalactic field -Bigm of 1CP 11 to 1CP 10 
G. The preceding amplification of the magnetic field in- 
side the dwarf galaxy by the combined action of stellar 
activity and differential rotation (i.e. the Galactic Dy- 
namo) is thereby restricted by the IGM magnetic field 
already present at the formation time of the galaxy. For 
an IGM magnetic field of -Bigm ~ 1CP 10 G, the Lorentz 
force may prevent the formation of a new generation of 
dwarf galaxies and subsequent star formation. As a con- 
sequence, the IGM magnetic field never grows signifi- 
cantly above 10 -10 G. Since dwarf galaxies are charac- 
teristic in the early phase of the evolution of the uni- 
verse, this Cosmic Dynamo may have been very efficient 
in magnetizing the IGM. However, besides the accretion 



of IGM material pr e vious ly enriched with magnetic fields, 
IDubois &: Tevs sicr (|2009f) also point out the importance 
of accretion of satellite galaxies for the evolution and am- 
plification of the magnetic field in galaxies at later times. 

Our simulations emphasize the Cosmic D ynamo sce- 
nario proposed by Du bois fc Tevssie r (2009). The effi- 
cient amplification of the magnetic field during the equal- 
mass-merger presented here clearly shows that interac- 
tions of galaxies should be taken into account in studies 
of the magnetic evolution of the universe. We would 
also expect an intergalactic medium which is not only 
enriched with heavy elements by stellar activity, but 
also magnetized on large scales by galaxy interactions. 
Our simulations may help to understand the observa- 
tionally well established facts that very young galaxies 
already exhibit magnetic field strengths comparable with 
nearby fully developed spiral galaxies and the rotation 
measure estimates of intergalactic magnetic fields (e.g. 
iBernet et al . 2008). 
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